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Drift velocity peak and negative differential mobility in high field transport in 
graphene nanoribbons explained by numerical simulations 
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We present numerical simulations of high field transport in both suspended and deposited armchair 
graphene nanoribbon (A-GNR) on Hf02 substrate. Drift velocity in suspended GNR does not 
saturate at high electric field (F), but rather decreases, showing a maximum for F ~ 10 kV/cm. 
Deposition on Hf02 strongly degrades the drift velocity by up to a factor of « 10 with respect to 
suspended GNRs in the low-field regime, whereas at high fields drift velocity approaches the intrinsic 
value expected in suspended GNRs. Even in the assumption of perfect edges, the obtained mobility 
is far behind what expected in two-dimensional graphene, and is further reduced by surface optical 
phonons. 
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Assessing the potential of Graphene Nanoribbons 
(GNRs) for future electronic applications requires full 
understanding of both quasi-equilibrium and far-from- 
equilibrium transport mechanisms 0, However ex- 
perimental low-field (LF) mobility in 1 nm-wide GNRs 
can be as low as 100 cm 2 /Vs and is limited by edge disor- 
der 0, Furthermore, we know that achieving ideally 
smooth edges is not enough: full-band (FB) modeling 
shows that LF mobility due to only acoustic (AC) phonon 
scattering at low fields is close to 500 cm 2 /Vs for 1 nm- 
wide GNRs [i[ . Most importantly, nanoscale transistors 
do not operate in the LF mobility limit. Therefore, simu- 
lation of far-from-the-equilibrium transport conditions is 
required to understand achievable device performance. 
Whereas at low field carrier scattering is mainly due 
to low-energy intravalley acoustic phonons [J], at high 
electric field scattering is dominated by optical phonon 
emission (EM), which becomes relevant when electrons 
gain enough energy to emit optical phonons. High-field 
steady-state transport can be simulated by solving the 
Boltzmann transport equation (BTE) through the single- 
particle Monte Carlo (MC) method. When dealing with 
ID systems, however, particular attention has to be paid 
in accurately describing the energy dispersion relation 
due to quantum lateral confinement: for each GNR sub- 
band there is a parabolic behavior close the subband min- 
imum, and then the typical graphene quasi-linear behav- 
ior already for relatively small wavevector values, corre- 
sponding to a velocity ps 8 x 10 5 m/s. Some authors 
have considered phonon confinement and multisubband 
transport, focusing on quantum wires with the effective 
mass approximation p|l{|. For GNRs, a multisubband 
MC approach has been followed in Q and BTE has been 
solved in a deterministic way at criogenic temperatures 
in Q. Bresciani et al. @ have used a 2D model, which 
is not fully adequate for sub- 10 nm GNRs, where size 
effects are indeed relevant. 

In this work, we adopt a steady-state single-particle 
full band MC approach accounting for carrier degener- 
acy [l(| , which has a significant effect for materials with 



a small density of states as graphene. Scattering rates are 
obtained within the Deformation Potential Approxima- 
tion (DPA) from phonon dispersions described by means 
of the fourth-nearest-neighbour force-constant approach 
(4NNFC) [ll| and a p z tight-binding Hamiltonian for 
the electronic structure. We consider in-plane longitu- 
dinal acoustic and optical (LA and LO), transversal op- 
tical (TO) and surface optical (SO) phonons. In each 
subband, the rates are computed on a 2000-point grid 
in the fc^-space (due to simmetry only longitudinal elec- 
tron wavevectors k x > have been taken into account), 
considering energy up to 1.5 eV above the bottom of the 
first subband and including up to 18 subbands: this en- 
sures accurate results even for strong longitudinal electric 
field F < 3 x 10 2 kV/m and for all the considered GNR 
widths (W < 10 nm). Due to Van Hove singularities in 
the ID rates, the self-scattering method [12( is inefficient, 
so that we have adopted the MC procedure described in 
Ref. [HJ and extended to quasi- ID systems [l3]. The 
whole story of an electron has a duration T = J^i ^-U 
(T < (10 — 100) fis depending on W and F), where 
At, = fj+i — U is the i-th time step and U is the i-th 
sampled time. At each ti, the wavevector k = (k x ,k yri ), 
where k yri is the transverse quantized electron wavevec- 
tor is computed and the average value of quantity 
X (either the drift velocity u(k) = 1/h ■ dE(k)/dk x or 
the energy EQt)) is evaluated according to Qjl]: (X) t = 

^ J T X [k x (t')] dt' = i J2 U ft! +1 x [M*')] dt ' ■ The elec " 
tronic temperature T e i under an applied homogeneous 
field F is computed by equating (v(k x ) 2 ) with the mean 
squared velocity at equilibrium [13j . For F = 0, T e i is 
equal to the lattice room temperature T/ at , whereas for 
F > 0, T e i > Ti at . Once obtained T e /, the distribution 
function is updated accordingly at each ti and final states 
for scattering are filled obeying to the Pauli exclusion 
principle (PEP). 

The computed average drift velocity (vd) limited by 
intrinsic phonons is plotted in Fig. QJi as a function of 
F for suspended sub-10 nm GNRs and for a carrier den- 
sity ri2D — 10 12 cm~ 2 . (vd) strongly varies within the 
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FIG. 1: (Color online) (a) Average drift velocity (yd) as a 
function of F (n2D = 10 12 cm -2 ) and (b) as a function of 
U2D for different W of suspended GNRs. In (a) MC results 
for zigzag CNTs [l4| with sub-10 nm circumference (2.7 nm 
and 8.7 nm) are also reported. In (b) F — 10 kV/cm. 

considered F interval and a maximum appears at high 
electric field (F = Fth)- For F < Fth, (vd) can be fitted 
by means of the Caughey-Thomas model, i.e. (vd(F)) — 
{HinfiF) /[l + {^inflF/vp)' 1 ] 1 /' 1 , where 7 pa (1.3-=- 2), /i ir[j o 
is the intrinsic LF mobility and w p is the peak velocity, 
which ranges from 2 x 10 5 to 4 x 10 5 m/s, depending on 
W (Table fl}. (w^) does not saturate with F, analogously 
to what has been observed in the case of zigzag CNTs 
[li], [Hj]. As shown in Fig. [TJi, (v d ) slightly decreases for 
increasing F [l4j]. However, in constrast to [IH, the neg- 
ative differential mobility (NDM) cannot be explained 
by the increased number of populated states in the sec- 
ond subband with smaller v(k). Indeed, as in [l5j|, we 
have verified for W = 10 nm that such effect is much 
more pronounced if we fictitiously limit transport to only 
one subband (red dashed line in Fig. QJi). The veloc- 



TABLE I: Caughey-Thomas parameters fii„fi (cm 2 /Vs), v p 
(m/s) and 7 for different W (nm), obtained through the fitting 
with MC data. 
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ity peak can be explained by the combined effect of the 
quasi-linear dispersion relation in GNRs and the strong 
increase of optical EM at high F, which increases the oc- 
currence of backscattering events therefore reducing the 
average velocity. If we follow the story of a single elec- 
tron, we can see that as F increases, instantaneous elec- 
tron velocity cannot increase beyond the limit imposed 
by the dispersion relation: any backscattering event will 
invert the instantaneous velocity sign and therefore re- 
duce the average velocity. Therefore, in the absence of 
optical phonon EM, the drift velocity (vd) would satu- 
rate to about 8 x 10 5 m/s. The onset of optical phonon 




Electric Field F (kV/cm) 

FIG. 2: (Color online) (a) Intrinsic mobility fii n and (Inset) 
average electron energy as a function of F. Results for GNRs 
with W = 5 nm 0] are also shown in (a). U2d = 10 12 cm -2 . 



EM makes (vj) peak at a fraction of that value and then 
decrease with F. We believe that the very same mech- 
anism explains also the NDM in zigzag CNTs, as well 
as in graphene [lt| , even if it has not been proposed be- 
fore uduM- 0ne 

can also see that the peak velocity v p 
increases with W. Indeed, if fku is the optical phonon 
energy (pa 160 meV for LA mode), the current J can be 
estimated as J = fojj [17], where vp is the Fermi 
velocity. Since v p — J jien^D) and the Fermi wavevector 
kF = nn2DW, we obtain v p = lu x W/ir oc W. 

The threshold field F t h strongly decreases with W, be- 
cause of the increased mean free path (Lk), which allows 
electrons to gain energy required to emit optical phonons 
at lower F. F t h can be roughly extimated by imposing 
the cutoff energy Hui for optical EM equal to the mean 
kinetic energy gained between two scattering events, i.e. 
qFth(Lk). Since by increasing W from 1 to 10 nm, (ik) 
increases from w 6 nm to 1 /im [4], F t h decreases from 
fa 2.7 x 10 2 kV/cm to « 1.6 kV/cm. We remark also that 
the obtained values for are in agreement with those 
found for zigzag CNTs with a circumference comparable 
with the considered W, both in linear and in non- linear 
regimes In Fig. [T]d, {vj) is plotted as a function of 
n 2 D for F = 10 kV/cm: (vd) does not depend on n 2 D 
even in the degenerate regime (n2D — 10 12 cm -2 ), where 
PEP limits up to the 50% of the scattering events. In 
Fig. [2 we show the intrinsic mobility jii n = (vg) j F as a 
function of F for n2D = 10 12 cm~ 2 . Since in the linear 
transport regime Vd oc F, fii n is constant for F < F t h 
and decreases above F t h as //j„ oc 1/F a with a > 1. In 
addition, the narrower the ribbons, the stronger the sup- 
pression due to lateral confinement. Results are in good 
agreement with a multi-subband model for W pa 5 nm 0] • 
In the inset of Fig. [21 the average kinetic electron energy 
(£k) in units of ksTi at is shown as a function of F for 
n 2 D — 10 12 cm -2 . While for low F pa 0.1 kV/cm elec- 
trons tend to remain near the first conduction subband 
edge and (i?k) ~ l/2fcsT/ at for the narrowest ribbons, 
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Electric Field F (kV/cm) 

FIG. 3: (Color online) (a) (vd) and (b) mobility flex cLS 3- 
function of F for GNR on Hf0 2 . n 2D = 10 12 cm -2 . 



for high field F > 10 2 kV/m, (£ k ) > fcT iat and higher 
energy states are occupied. Note also that (Ey) increases 
with W, since subbands become closer, allowing electrons 
to populate higher subbands. 

In Fig. [3^, and|3}3 the average drift velocity (vd) and 
the mobility /j, ex = (vd)/F are shown, respectively, as a 
function of F for GNRs deposited on Hf0 2 including 
the effect of SO phonons (n 2 D — 10 12 cm' 2 ), through the 
first SO(l) and second SO(2) modes {frw so = 12.4 meV 
for SO(l) mode). With respect to the suspended GNR 
case, (vd), as well as mobility in GNR on Hf02 are almost 
one order of magnitude smaller in the LF regime, while 
similar values are obtained for high field, as already noted 
for graphene on IHO2 [Hj]. As in graphene [l6j], deposi- 
tion on HTO2 leads to an extension of the linear region to 
fields up to a factor 10 larger than those corresponding 
to suspended GNRs. For narrow ribbons, (vd) does not 
saturate even for F — 3 x 10 2 kV/cm. 

In order to understand the different behaviour in sus- 
pended and deposited GNR on HfC>2, in Fig. 0k we show 
the distribution function G as a function of (Ey — Eqi), 
where Eci is the first subband edge, for F = 1 and 
10 kV/cm (W = 4.86 nm) for GNRs, both suspended and 
deposited on Hf02. At low fields, G decreases rapidly 
with energy in both cases, showing, unlike graphene, 
sharp peaks due to intersubband scattering. At high 
fields instead, electrons are excited up to energies close or 
above 1 eV, increasing T e i (inset of Fig. [2]). As expected, 
deposition on Hf02 leads to a shorter high energy tail 
compared to that in intrinsic GNR (Fig. 0Ji) due to the 
introduction of an additional channel for energy and mo- 
mentum relaxation which increases total scattering rate 
and pushes electrons down to lower average electron en- 
ergy. As for graphene on HTO2 [III, the rough absence 
of degradation of (vd) at high F can be explained by 
the reduced population in the high energy tail, i.e. by a 
smaller amount of electrons in the nonlinear dispersion 
energy region where the band velocity is smaller, which 
more than counterbalances the decrease of (vd) resulting 
from the increased total scattering rate when deposit- 
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FIG. 4: (Color online) (a) Energy distribution function for 
F = 1 and 10 kV/cm for both suspended GNR (continuous 
lines) and GNR deposited on Hf02 (dashed lines). Fraction of 
different scattering events as a function of F for (b) suspended 
GNR (LA, LO and TO phonons) and (c) GNR deposited on 
Hf0 2 (LA, LO, TO and SO phonons). W = 4.86 nm and 
n 2D = 10 12 cm -2 . 



ing GNR on substrate. In order to explain the increase 
of the F interval where (vd) shows a linear behavior in 
the GNR on Hf02 case, in Fig. |4]d and |4j: we show the 
relative ratio of scattering events for the different mech- 
anisms as a function of F for 5 nm-wide suspended GNR 
and GNR on HTO2, respectively. For intrinsic GNR and 
F < 1 kV/cm, the main scattering events involve ab- 
sorption (ABS) and EM of LA phonons (Fig. gp). For 
F > 10 kV/cm, optical phonon EM (of both LO and TO 
phonons) becomes predominant, increasing up to 30% of 
the total number of events. For GNR on Hf02 instead, 
scattering involving SO(l) phonons happens very often 
already for F = 1 kV/cm whereas rates of LA, LO and 
TO phonons are limited to few percents even for high 
fields, as can be seen in Fig. |4j:. In particular, the ob- 
served large absorption rate of SO(l) phonons (Fig. Eh) , 
which is associated to the high Bose-Einstein occupation 
factor, appears to counterbalance SO(l) emission even at 
high F, and is the responsible of the extension of the lin- 
ear region up to fields of 10 kV/cm or above, depending 
on W. 

In conclusion, we have performed a FB investigation 
of the dependence of drift velocity and mobility on the 
electric fields in GNRs. Suspended GNRs exhibit a drift 
velocity peak and then a NDM for large electric field, 
as also observed in zigzag CNTs. This property is due 
to the combined effect of quasi-linear dispersion relation 
and the emission of optical phonons. In particular, the 
maxima occur for a threshold field F t h oc 1/ (Ly) oc 1/W. 
Depositing GNR on HTO2 substrate strongly degrades 
(vd) at low field by a factor w 10, whereas at high fields no 
degradation is observed [llj]. Furthermore, in deposited 
GNRs velocity saturation and peak are shifted at higher 
F, due to the compensation of SO absorption and EM 
mechanisms. 
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